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Summary 


When an energetic electron is incident on a solid surface, a considerable 

rnmoonents such as solar cells In low earth orbit. The purpose of this report is to 
pSe a convenient scheme for including realistic SEE behavior in numerical 
simulations of plasma-suuriace interactions. 

Previous Monte Carlo simulations provide a data base for properties of 
/■'nnriarv plectron emission (SEE) from insulators and metals. Incident primary 
S5 ensure considered at energies up to 1 200 eV. The behavior of 
electrons is characterized by (1) yield vs. pnmaiy energy E p . (2) distribution vs. 

secondary energy E s and (3) distribution vs. angle of emission 9. 

Fnr nrimarv eneroies above 50 eV the SEE yield curve can be conveniently 
nnrflmptprized bv a Haffner formula. Special attention is paid in this report to the low 
^"nge^ E ^p to 50 eV where th P e number and energy of secondary electrons 

is limited by the finite band gap of the insulator. There is also a considerable 
DTObability for elastic backscattering or reflection of the pnmary electron in t 
energy range, and separate yield curves are given for reflected primary electrons. 

The energy distribution of secondary electrons is described by an empirical 

ir & - n 

slightly more peaked in the forward direction than the customary cos 0 dist 
The results of the MC simulations-yield and distribution curves-are 

solid surface interactions in low earth orbit. 
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I. Introduction to Secondary Electron Emission 


A. Description of SEE Processes and Behavior 


When electrically charged particles with sufficient kinetic energy are incident on 
the surface of a solid material, the material emits electrons. This paper will only con- 
sider electrons as the incident charged particle. In keeping with the current literature 
these incident electrons will be called “primary electrons,” and the electrons from the 
material that are excited by the primary electron will be called “secondary electrons.” 
The secondary electron yield 8 is the ratio of external secondary electron current to pri- 
mary electron current. The yield may also be regarded as the average number of ex- 
ternal secondary electrons per primary electron. 

Different types of materials exhibit characteristically different yields. The basic 
source of the difference in yield between the types of materials is the electron energy 
band structure. A diagram showing the energy bands for a metal and for an insulator 
is shown below. 


Metals Insulators 



Fig. 1 Energy level diagrams for metals and for insulators. 
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The single headed arrows represent the amount of energy an internal electron 
may acquire by interacting with the primary electron or other internal secondary elec- 
trons. As indicated in the diagram, if the internal secondary electron has sufficient en- 
ergy to climb to the vacuum level, it may be able to leave the material and become an 
external secondary electron. It is important to keep in mind, though, that the internal 
secondary electron must also be traveling toward the surface of the material if it is to 
become an external secondary. On its way to the surface, an internal secondary may 
interact with many other internal electrons. With each interaction the secondary elec- 
tron will lose energy and/or suffer a change in direction. It is the number of these inter- 
actions and the distance to the vacuum level that ultimately determines whether or not 
an electron will become an external secondary electron. 

In metals the work function <D is the energy that all electrons must be able to 
overcome to become external secondary electrons. Just below the work function is the 
top of the conduction band E F . The conduction band is partially filled, with occupied 
levels just below E F and unoccupied levels immediately above Ef. The core represents 
those energy levels that are already completely filled. These are the only energy lev- 
els in which an electron can exist. Therefore, if a primary electron penetrates the sur- 
face of a metal with sufficient kinetic energy, it will be in the conduction band since it is 
not able to enter the core by the Pauli exclusion principle. 

Since the conduction band is partially filled, both the primary electron and the 
internal secondary electrons will often suffer inelastic collisions in which only a small 
amount of energy (a few eV) is lost to conduction band electrons. These frequent 
small-energy-loss collisions cannot lead to external secondary electrons since the en- 
ergy transferred between the colliding electrons is less than O. Therefore much of the 
initially available energy is “wasted” and the yield is small. 

The situation is quite different for insulators. Insulators have an electron affinity 
X that all electrons must overcome to become secondary electrons, x ' s the energy 
difference between the surface of the material and the bottom of the empty conduction 
band. Below the conduction band is the valence band, which is completely filled for 
good insulators. The core levels are at still lower energies. Between the bottom of the 
conduction band and the top of the valence band is the band gap E g ; no electrons are 
allowed in this energy region. In insulators the width of the band gap is the major de- 
terminant of the secondary electron yield. 

As illustrated schematically by the diagram on page 4, a good insulator is char- 
acterized by a large band gap. In this case most excitation processes raise an internal 
electron above the vacuum level, allowing for possible secondary electron emission. 
This is the basic reason why the maximum yield for insulators is much larger than for 
metals. Semi-conductors have a small band gap. For such materials the situation is 
similar to that of metals; hence the maximum yield for semiconductors is between 
those for metals and insulators. For very large band gap insulators the maximum yield 
tends to decrease somewhat, perhaps because the primary electron can penetrate to 
greater depths before exciting the internal electrons. 
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Small Band Gap 
(semiconductor) 



Large Band Gap 
(insulator) 



Fig. 2 Energy level diagrams for a semiconductor and for a large band gap insulator. 


The diagram below shows the maximum yield versus the band gap for several 
real solids (open circles), and the results of G. A. Rezvani’s Monte Carlo simulations 
(solid circles). The materials are identified in the table on the following page. 
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Fig. 3 Maximum SEE yield versus band gap E g (eV) for a number of solids. 
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No. 

material 

E g (eV) 

X (eV) 


E m (eV) 

i 

P bS 

0.4 

4.7 

1.25 

500 

2 

Sb 7 Te 3 

1.1 

4.34 ' 

1.2 

700 

3 

GeS 

1.8 

4.0 

1.0 

400 

4 

PbO 

2.6 

3.6 

1.8 

600 

5 

Te0 2 

3.0 

3.4 

1.7 

600 

6 

S b^Oz 

4.2 

2.8 

1.6 

500 
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Nal 

5.8 

1.5 

8.35 

1390 

8 

KJ 

6.2 

. 1.1 

11.68 

1430 

9 

C si 

6.3 

0.1 

17.23 

2150 

10 

CsBr 

7.0 

0.2 

18.61 

2340 

11 

CsCl 

7.5 

0.3 

14.68 

1470 

12 

NaBr 

7.7 

- 

10.34 

1430 

13 

KBt 

7.8 

0.9 

8.35 

630 

14 

RbCl 

8.2 

0.51 

10.62 

1310 

15 

LiBr 

8.5 

1.0 

6.22 

580 

16 

TUI 1 

8.5 

0.6 

10.14 

1000 

17 

NaCl 

8.6 

0.65 

11.39 

1300 







18 

LiCl 

10.0 

0.54 

8.14 

1000 

19 

NaF 

10.5 

- 

12.38 

1290 

20 

KF 

10.9 

0.1 

10.7 

1100 

21 

LiF 

12.6 

1.0 

6.42 

540 


Table 1 Band gap, electron affinity, maximum yield, and energy of primary 
electron at maximum yield for several insulators. 


The torturous nature of the SEE process is illustrated by the diagram on page 6 
from one of Rezvani’s simulations. The solid line shows the path, projected onto the 
yz- plane, of a 400 eV primary electron incident normally on solid argon. The dotted 
lines show the paths of the internal secondary and tertiary electrons. (A tertiary elec- 
tron is one excited by a higher energy internal secondary electron, etc.) According to 
Rezvani, the primary electron suffered 200 interactions (the maximum allowed by his 
simulation) and generated 22 internal secondary electrons and 7 tertiary electrons. In 
this simulation run 12 of the internal secondary electrons reached the surface and 
managed to escape. 



•sauji U9>)0jq hum 

S8UEDUOO0S 0MI *0 S8U0P8fe4 eqi PUE ‘8UIJ P! |0S e qiiM UMOMS SI uojpeie 
XjBiuud 8Mi *0 Ai 0 P 8 fEj; 8M1 -BUBid-zA ei» opo AjopefBn eqi P uoipe wd 

am SMOMS 0jn6|i 8M1 uo6je ui. uojp0|0 AjBiuud OOfr peeJiAuuJBj Bid 






7 


B. Characterization of External Secondary Electrons 


The most important characteristics of secondary electron emission are the sec- 
ondary electron yield versus the primary electron energy, the energy distribution of the 
external secondary electrons, and the angular distribution of the secondaries. 


Most of the research has focused on the first characteristic, secondary electron 
yield as a function of primary energy. It has been shown that for all bulk solids the sec- 
ondary electron yield will start from zero at E p = 0 eV, rise rapidly to a maximum yield 
5 at E m , then slowly decrease with E p > E m . If the secondary electron yield for some 
materiaT were plotted versus energy of the primary energy, the result would be a varia- 
tion of the “universal” yield curve. Such a curve is shown below. 



Fig. 5 Universal SEE yield curve. 

Typical values of the maximum yield for various materials are shown in the fol- 
lowing table. Notice that the yield is greater than one, often by a large amount, for 
semiconductors and insulators. 

Material Maximum Yield 

Metals less than 1 

Semiconductors 1 * 2 

Polymers 2-3 


Insulators 


6-12 



It is convenient to define E, and E 2 as the energies of the primary electron for 
which the yield is 1.0. These energies are significant for plasma-insulator interactions; 
for E 1 < E p < E 2 the yield will be greater than 1 .0, and the insulator will charge posi- 
tively. Likewise, for Ep < E 1 and E p > E 2 the yield will be less than 1 .0 and the ma- 
terial will charge negatively. For insulators E : is about 1 6 eV, is around 200 - 400 
eV, and E 2 is typically greater than 1000 eV. 

The next important characteristic of secondary electron emission is the energy 
distribution of the secondaries. The major features of this distribution are illustrated on 
the next page for a typical primary energy E p = 150 eV. The maximum S at 10 eV is 
nearly independent of E p , and is therefore representative of “true” secondaries. These 
are electrons that were initially present as bound electrons in the solid, were excited 
by the primary electron, and finally managed to escape from the surface to become ex- 
ternal secondary electrons. 

The maximum B is located at the primary energy E p , and can only be from pri- 
mary electrons that have undergone elastic back-scattering. The shape of the small 
maximum U is characteristic of the material and originates from inelastically scattered 
primaries, i.e. primary electrons that suffer only one or two inelastic collisions before 
being scattered back out of the solid. The region between 50 eV and the maximum U 
originates from a mixture of true secondaries and inelastically scattered primaries. 

The two types of secondary electrons in this region are not distinguishable by means 
of external measurements. 

The results of secondary electron emission described in Part II of this report are 
based on Monte Carlo simulations of the SEE processes. To simulate the elastic scat- 
tering, a general overall form of the elastic scattering cross section typical of low atom- 
ic mass elements (12 to 40 amu) was used. Since detailed features of elastic scatter- 
ing for specific elements were not included, the simulations were not able to reproduce 
the back-scattering peaks at U and B for primary energies E p > 50 eV. However, for 
primary energies 0 < E p < 50 eV the back-scattering peak at B was obtained. In all 
cases the simulations give good results for the main secondary peak at S, which is 
more important for the purpose of modeling the plasma-insulator interactions. 

The final important characteristic of secondary electron emission is the angular 
distribution of the external secondaries. Because of the large number of collisions suf- 
fered by the primary and secondary electrons inside the solid (tens or hundreds in typi- 
cal cases), the angular distribution of excited electrons near the surface should be es- 
sentially isotropic. But since the excited electron must be traveling toward the surface 
to be emitted, the angular distribution of external secondaries should approximate a 
cosine weighted distribution, i.e., weighted in the direction of the outward normal to the 
surface. Rezvani found that this distribution is essentially independent of the energy 
of the primary electron and of the angle of incidence of the primary electron. A straight 
forward cosine distribution is generally accepted for experimental results. However, 
the Monte Carlo simulations give a distribution peaked slightly more in the forward 
direction, as will be presented in Part II. 
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Energy (eV) 


Fig. 6 Schematic energy distribution of secondary electrons. This diagram is not 
drawn to any particular scale and does not represent any particular material. 
However, the broad features of the diagram are true of almost all materials. 


In these simulations, Rezvani observed that after a primary electron, or excited 
internal electron, suffered two or three interactions with other electrons it would “forget” 
about its initial direction of travel. Furthermore, most internal secondaries are excited 
with modest energies of 5 to 50 eV, and after excitation become “ignorant” of the ener- 
gy of the primary electron. This behavior appears vividly in the jagged nature of the 
tracks for both primaries and secondaries (Fig. 4). Consequently, both the energy dis- 
tribution and the angular distribution of secondary electrons are nearly independent of 
the primary energy and of the angle of incidence of the primary. 

This independent behavior on the part of the secondary electrons is observed in 
all bulk solid materials. For very thin films, however, more detailed features of the 
paths would have to be taken into consideration. From Fig. 4 the mean penetration 

O 

depth of the primaries is about 50 4, and most secondaries do not wander to depths 

v a 

below 200 A. Therefore, “thin" means only a few hundred Angstrom units. We shall 
concern ourselves with bulk materials in the solid state, and refer the interested reader 
to the literature for thin films. 

In the next part of this paper we will develop a method of simulating secondary 
electron emission for use in numerical simulations of plasma-insulator interactions. 
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II Method of Including SEE in Numerical 

Simulations of Plasma-Insulator Interactions 


When an insulating solid is situated in a plasma environment, the resulting in- 
teractions affect The properties of both the insulator and the plasma. The insulator will 
acauire a surface charge which may be negative or positive depending on the poten 
Erf the "or She properties of the plasma. In addition the density r and vetoc 
tv of the particles in the plasma will be affected out to a distance of several Debye 
enothJ aww from the insulating surface. One way of studying the behavior and 
Dhvsics of these interactions is to create a computer simulation to follow the motion of 
the charged particles in the volume of the plasma near the surface In any such pro 
cram an^mportant feature of the plasma-insulator interactions is the production of 
secondary electrons when an electron from the plasma impacts the surface of the 
fnsulator The pur p 0Se of this report is to describe a method for incorporating he 
SEE response of an insulator (or metal) surface into a computer simulation of the 
plasma-insulator interactions. 


The situation we wish to describe goes as follows. An electron in the P la ^ma. 
moving under the influence of the electric fields produced by all the charges in the sy 
tern imoacts at a certain time somewhere on the surface of the insulator. From 
locitv components of this primary electron, we can determine its kinetic energy an 
angle ofSence onto the surface. At this point in the simulation, we need to deter- 
mine how many secondary electrons are emitted. If there are some, we wiH need to 
determine the kinetic energy and direction of motion of each of the secondaries. In th 
following sections of this part we describe the distributions and procedures tha 
used for this purpose. 


A. Total Yield vs. Primary Electron Energy 


The secondary electron yield 6 depends on the kinetic energy and angle of me 
Honrp nHhe orimarv electron and is also influenced considerably by the properties o 
fhTmaterianhaUoms fhTsurface. Typical behavior .or mode, quaru and alornmum 
is shown on the figures on page 14 . At zero primary energy “ ' s r “™JS^ lh , 
even for primary energies of a few eV there are a considerable number of elastically 
reflected primaries; see section D.) As the primary energy increases the yie'd 'n- 
ae e ase e s d r P a r idiy, simply beoause more energy 

The yield reaches a maximum at some energy E m , and then decreases siowiy ior y 
er nrimarv energies. This is because higher energy primaries penetrate more deeply 
into the material which implies that the excited electrons have further to travel before 

escaping from the surface. 
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A convenient fitting formula for 5 0 , the yield at normal incidence (9 = 0), is pro- 
vided by the Haffner formula. 


5 0 (E P ) = C (e-^ - e- bE p) 


The parameters a, b, and C are empirically fitted to experimental or simulation data for 
a particular material. To relate 5 m and E m to these parameters, we start by setting the 

derivative of 5 equal to 0: 


Q _ 95 0 (Ep2 = c L ae -aE p + be- bE p) 
oEo ' 


( 2 . 2 ) 


ae-aEp = be _bE p 


(2.3) 


ie 

>e 



= lp g (keV) 
(b - a) 


(2.4) 


This is the energy of the primary electron for which the yield 8 is a maximum. 
Assuming the maximum yield 8 m is known from the data, the parameter C can be 

found by inserting E m into the Haffner formula: 


ci- 

of 


jh- 

y 



=> c = 



(2.5) 


( 2 . 6 ) 


For the plasma-insulator simulations we should like to use quartz as a prototype 
for a good insulator. But the complexity of the composition and crystal structure of 
quartz prevent us from a direct simulation of the secondary electron emission, and the 
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experimental data are old and questionable. Therefore to find parameters for quartz 
we look at our simulation results for KCI, which has a band gap of 8.5 eV that is 
comparable to 8.0 eV for quartz. Similarly, the electron affinity for quartz is 1.0 eV, 
whereas the electron affinity of KCI is 0.6 eV. We have performed extensive simula- 
tions with model KCI for X = 0.0 eV and X = 2.0 eV. Since the electron affinity for 
quartz is halfway between these values, we propose a simple linear interpolation be- 
tween the results for the two models of KCI, as illustrated below. 



Ep (keV) 


Fig. 7 Haffner yield curves for several model solids. 

(a) Model KCL with X = 0.0 eV 

(b) Model quartz with X = 1 .0 eV 

(c) Model KCL with X = 2.0 eV 


The secondary electron yield also increases with increasing angle of incidence 
of the primary electron, where 0 is measured from the normal to the surface. As 0 
increases, the depth to which the primary electron penetrates decreases; thus the ex- 
cited internal electrons have less distance to travel to the surface, and the secondary 
yield becomes larger. An empirical formula to represent this effect is: 
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6(0) = S 0 e 


c(l-COS 0) 


(2.7) 


Some references suggest a value of c = 2 for polymers, but this choice would make 
S(90°V5 = e 2 = 7.4. From a few Monte Carlo simulations of SEE at fairly large 
angles (8 = 60° to 80°), our experience indicates that this factor is too large. We prefer 
to use c = 1 , which gives 5(90°)/5 o • e = 2.7 instead. 

On page 14 are plots of the Haffner formula for quartz and a'uminum w-ith 
angle 6 varying from 90° for the greatest yield down to 0 in increments 
values of the parameters we used are shown below. 


Material 

Quartz 


Aluminum 


Parameter 

a = 1 ,02/keV 
b = 5.47/keV 
C = 14.02 

a = 0.5/keV 
b = 8.3/keV 
C = 1.23 


Em 


377 eV 


360 eV 


Monte Carlo simulations, as shown in Fig. 7. 

in a numerical simulation that would involve 

ables for the Haffner formula are calculated as follows. 


ze 


x- 

y 


E p 

0 P 


nyTW+V? + V? 

2 x 1.6021 89xl0- 19 


= arctan 




V 


V, 


eV 

radians 


( 2 . 8 ) 

(2-9) 



SECONDARY ELECTRON YIELD 


V.L e'UU.U 400.0 DOO.D 800.0 1000.0 


At the maximum yield is 21 60 

At 6 if the maximum yield is 12 79 

At 3Cf Lhe manmum yield ls 6 8" 

At 0 (f the maiimuEr. yield is 7 76 


ENERGY OF THE PRIMARY ELECTRON (eV) 


1200.0 



At 90* the maximum yield is 2 61 
At 60'the maximum yield is 1 56 
At 30°lhe maximum yield is 1 10 
At 00° the maximum yield is 0 96 


ENERGY OF THE PRIMARY ELECTRON (eV) 


Fig. 8 Haffner yield curves for various angles of incidence 
of the primary electron. 
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The yield 5 from the Haffner formula will be a fractional number. This number is 
a ratio of primary electron current to seconda^ electron current. However, ,n a burner- 

to be the average number of secondary electrons emitted from the matenal per mci 
dent primary electron. 

The method we use to accomplish this result te 9 ins ‘^^^ihe 

VIP1 nMAX which calculates the Haffner formula as in equation (2.7) and returns t 
S D y" men add a random fraction to the yield and truncate the sum to an 

inte ^ er ‘ y = RAN(seed) 


N 


pos 


= INT(yield + y) 


N . More on this later in section C.) 

Sffhefi " be less than 2. Repeating the process for a large number of secon- 
daries, an average value of 1 .7 will be obtained for the yield. 


B. Energy Distribution of Secondary Electrons 


Each secondary electron that is produced in a numerical 
assigned an energy and a direction. At present our concentration win be on assigning 
an energy to each secondary electron with a suitable distribution funct 

We developed our distribution function based on Rezvani’s simulation data tar J 00 
n rimary electrons with E p = 50.0 eV at normal incidence on solid argon. The his 

togram of these results is shown in the figure on page 18. Since the [ er ° 

secondaries climbs rapidly at low secondary electron energies and falls off rap y 
. i- a § _i : i i. i.i * i+inn fi inrtinn of the TOrm. 


COflUaMt;i> OMllluo lapruiy at w , 

increases past 5.0 eV, we decided to try a distribution function of the form. 
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f(EJ = A VE" e'® 5 *”’ 2 

=5- dP(E s ) = f(E s ) dE s 


1 i 


Now we can establish some properties of this distribution function. First we nor 
malize the integrated probability to unity. 


1 = P(E S ) = I dP(E s ) = I f(E s ) dE s 

0 JO ( 2 . 12 ) 


/• oo 


= A 


Va e-^ 2 dE s 


o 


( 2 . 13 ; 


Making the substitution u = (E s /E 0 ) 2 produces an integral that can be expressed in 
terms of the gamma function r(u). 


1 = 


_ AEf 


u 3/4 • 1 e* u du 


JO 

( 2 . 14 ) 

\ A Ef no.75) 

( 2 . 15 ) 

A = 1.632 E 0 3/2 

( 2 . 16 ) 
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To find the maximum of the distribution 
equal to zero: 


function, we set the derivative with respect to Ej 


0 = 


9f(E s ) 

3E< 


E ^ 2 _ X jc - 1/2 

Eo 4 


(2.17) 


Eo 


>max 


( 2 . 18 ) 


The average energy is found from 


M OO 


p OO 


E a v — 


E s f(E s )dE s = A 


o 


E 3/2 e -(Es/Eo) 2 dEs 


0 


(2.19) 


3) \ 

With the usual substitution u = (E s /E 0 ) 2 the result is 


E 


av 


ra.75)^ 

-E'O 

r(0.75) 


0.9064 - 

1.2255 



0.7396 E 0 

( 2 . 20 ) 


A least squares fit of this distribution to Rezvani’s data gives E 0 10.28 eV 
! 4 ) Rezvani’s data and the fitted distribution function are shown on the next page. 


15) 


16) 
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0.0 5.0 10.0 15.0 20.0 25.0 30.0 


ENERGY OF THE SECONDARY ELECTRON (eV) 

Fig. 9 Energy distribution of secondary electrons. The histogram gives the results 
of the Monte Carlo simulations for E p = 50 eV in argon. The curve shows 
the best fit of the empirical formula of equation (2.10). 


Notice that this distribution is ngla Maxwell-Boltzman energy distribution. In 
fact, early attempts to fit the histogram with a MB distribution gave poor results. If the 
peak of the MB distribution is chosen at about 5 eV to fit the histogram, then the high 
energy end of the MB distribution is much too large for the results of the distribution as 
displayed on the histogram. Physically, there is no reason to expect the secondary 
electrons to fit a MB distribution, which corresponds to thermal equilibrium. Although 
the secondary electrons make a number of collisions with atoms in the solid before es- 
caping (the number of collisions can be anywhere between one and several hundred), 
they do not by any means come into thermal equilibrium with the solid or with each 
other. Consider, for example, the average energy of the secondary electrons at 7 5 eV' 
this is much greater than kT room = 0.025 eV. 

It will be seen that most of the secondary electrons have modest energies, in the 
range 5 to 1 5 eV. This is significant for the plasma-insulator simulations in two ways. 
First, the secondary electrons have enough energy to wander away from the surface 
for a significant distance (several centimeters) which allows them to drift in the electric 
field present near the surface. But when these secondary electrons fall back onto the 
surface, they produce a yield less than one, as given by the Haffner formula eq. (2.1). 
However, as will be discussed in section D of this part, at the lower energies there is a 
significant probability that the impacting electron will be elastically reflected. These 
effects can produce a significant current sheath just above the surface of the insulator, 
which can be considered to be a surface current. Without secondary electron emis- 
sion this effective surface current would not be present. 
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Once the distribution function f(E s ) is known, there are two more steps in order 
to select energies in accordance with the distribution function. First, we need to obtain 
the integrated probability from 0 up to each E s 


P(E S ) 


rE s 

f(E s ) dE' s 

Jo 


( 2 . 21 ) 


The values of P(E S ) increase monotonically from zero at E =0 up to unity at E - 

Actually, we cut off the distribution at E s = 30 eV because the values of f(E s ) are 
negligible for E s > 30 eV. Second, we need to select a random number 0 < r < . 
With this number we find the value of E s for which 


P(E S ) = r (2.22) 


The solution of this equation can be written symbolically in terms of an inverse func- 
tional relationship. 



(2.23) 


Sometimes it is possible to carry out both steps of this process analytically, and 
then the last equation can be written as a specific algorithm to give E s directly from r 
This state of affairs occurs, for instance, for a two dimensional MB distribution (but not 
for one or three dimensions!). In our case the integrated probability would involve an 
incomplete gamma function; we have yet to come across the inverse ol a Junc- 
tion Our alternative is to integrate the distribution function numencally and^akea 
table of P(E ) vs. E s . Then to economically accomplish the solution of equation (2.22), 

we use interpolation to construct the inverse table, namely E s vs. P(E S ), at uniform 
intervals of P(E S ). We can then quickly solve equation (2.23) by table look-up and in- 
terDolation Furthermore, since both tables have as their first column evenly spaced 
entries we can use a one line algorithm to find the entries we want without searching 
through the tables. For example, say we have a value for P(E.) and want to find the 
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corresponding E s . Assuming there are 100 evenly spaced entries for P(E S ), which 
ranges from 0 to 1 , we find the index I as 

I = INT(100*r) + 1 (2.24 

Then the desired value of E s lies between E S (I) and E S (I+1), and we can quickly 
find E s by linear interpolation. 

The procedure we have developed in this section is based on the assumption 
that the energy distribution for E s of the secondary electrons is independent of the 
energy E p of the primary electron. For E s > 50 eV this assumption is reasonably 
valid. Once an internal secondary electron is excited, it has no further knowledge of 
E p , and most internal excitations occur in the range from 0 to 40 eV. We have also 
checked this assumption with the results of the Monte Carlo simulations of SEE for 
argon with E p = 50 eV and E p = 400 eV; there is no discernible difference between 
the two histograms for the secondary electron energy distribution (although there is a 
sizable difference in yield for the two cases). When E p < 50 eV, however, less ener- 
gy is available for formation of secondary electrons, and additional considerations 
must be take into account. These considerations will be considered in the next sec- 
tion. There is also a much larger chance for reflection of the primary electron as 
discussed in section D. 
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C. Considerations on Yield and Energy of Secondary 
Electrons at Low Primary Electron Energies 

We can now look in detail at the energy considerations involved as secondaries 
are liberated from bulk solid material. Consider the energy level diagram below. 
When a primary electron strikes an insulator with some energy E p , it may drop at most 
to the bottom of the empty conduction band since all of the states in the valence band 


Yacuum 
level — 



The bottom of the conduction bend 
is ell the further the primery elctron 
can penetrate. At this point the primery 
electron h8s lost an amount of energy 

AE max = Ep + X 


Conduction 
Band 

Electron affinity. 

Band gap 
V energy. Eg 

< Region that electrons 
cannot stay in 



Valence 

Band 




Core 


Fig. 10 Energy level diagram for an insulator. 
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are filled. The most energy that can be delivered to the material by the primary elec- 
tron, then, is the sum of its kinetic energy and the electron affinity X: 


AE max ~ + X 


(2.25) 


Any internal electron that becomes sufficiently excited because of the energy 
liberated by the primary electron may become a secondary electron. Each secondary 
electron will come from either the valence band or the core. However, the deeper the 
excited internal electron is, the more energy it will require to get to the surface of the 
material and become a secondary electron. The most loosely bound electrons in the 
solid are at the top of the valence band. These will require a minimum energy: 



(2.26) 


just to reach the vacuum level. E b is the surface barrier. Any excited internal electron 
having less energy than this will have no chance of becoming an external secondary 
electron. 

On the other hand, if an excited internal electron is able to make it to the surface to be- 
come a secondary electron, the greatest energy it can possibly have is: 


E Smax ~ AE max - Eb 

= (Ep +x) - (E g + X) 
Es max = Ep - Eg 


(2.27) 


Now it is necessary to insist that Es max > 0 in order to have any production of 
secondary electrons. Hence there is a minimum value of E p> which is obtained by 
setting Es max = 0. 


E 


= Ep . = Eg 

Emin fc 


(2.28) 
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It a primary electron strikes the surface with E p < E 0 , then the yield of secundary elec- 
trons must be zero. There is still the possibility, however, that the primary electron can 
be elastically scattered back out of the solid. The treatment of reflected electrons will 
be discussed in section D. 

We can now consider the procedure for assigning an energy to the first sec- 
ondary electron when E p > E 0 . We calculate 5(E p ) from the Haffner formula, equation 
(2.7), and then obtain the integer value N pos by the algorithm on page 1 5. The maxi- 
mum energy available for the secondary electron is given by equation (2.27). 

E s , = Ep -E g 

^max 1 k & 


Therefore we must truncate the energy distribution function of section B at Es^xi. The 
total probability that the secondary electron will have an energy between 0 and Es max i 
is 1 00. In order to obtain this result, we must renormalize the distribution function be- 
tween 0 and Es max as follows. We have a table of integrated probability for the com- 
plete distribution function according to 


r E s 

P(E S ) = p(x) dx 

Jo 


(2.29) 


By table look-up from this table we obtain 


P 


max 


PO^Smax l). 


(2.30) 


Then we define the truncated distribution by 


p'(E s ) 


p(E s ) 

P 

1 max 


and the same relation follows for the integrated probability. 

P'(E S ) = 

*max 


(2.31) 
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It will be seen that the integrated probability up to and including Es max 1 is unity, as it 
should be. 

In order to select an energy E s1 for the first secondary electron, we choose a 

random number r such that 0 < r < 1 . We want to set r equal to the integrated probabili- 
ty for some E s . 


r = P'(E S ) 


P(E S ) 

p 

1 max 


(2.32) 


So we want to find E s such that 

P(E S ) = T P m ax (2.33) 

We can find the desired value of E s by table look-up and interpolation in the inverse 
table for E S (P), just as in the preceding section. The advantage of this approach is that 
we do not need to actually construct a new table of P^Eg) for each value of Es max . 
We simply look for a different value of P(E S ), according to equation (2.33). 


When 5(Ep,0) is large enough to give N^g > 2, there is the possibility of produc- 
ing a second ( and perhaps a third, etc.) secondary electron. But there is only a finite 
amount of energy to produce secondaries and enable them to overcome the surface 
barrier. So we need some way of ensuring conservation of energy. We can do so with 
the following equations. 


AEmax 2 ^Emax " E s j ~ E-b 


(2.34) 


= Ep + % - E sl - (E g + %) 




Es max 2 “ AE max 2 (E g + %) 

= Ep - E sl - 2 E g - x (2.35 

^ s max 2 - 0- we d 0 n °t have enough energy to form another secondary electron, and 
we do not allow for reflection of the primary electron in this case, so we quit. But if 
Es max2 > 0. we proceed to select E s2 by the method of the preceding paragraph. 
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When N pos > 3 , we continue to look for the possibility of producing more secon- 
daries The above process continues until either the number of secondaries produced 
reaches N pos , or Es max j < 0. The actual number of secondary electrons produced 

ig called Ng^. Of course in all circumstances, — Npog’ 

We can use the same procedure when dealing with secondary electron emis- 
sion from metals. However, metals do not have a band gap nor an electron affinity. It 
is the work function <J> that comprises the surface barrier in metals. To correctly treat 
metals, replace x ' n equation (2.25) and E g + x * n equation (2.26) with <t>. Then 

equation (2.27) becomes 

Es max = AEmax “ E b 

= (Ep + <t>) - <|> = E p ( 2 . 36 ) 


Thus it is possible to form secondaries in metals for any E p > 0; there is no threshold 
energy. The reason for this difference in behavior from insulators is due to the differ- 
ence in band structures. In a metal there are unoccupied electron energy levels imme- 
diately above the Fermi level, so it is possible for the primary electron to drop right 
down to the Fermi level. In insulators, this behavior is not possible because of the 

finite band gap. 

For most metals the maximum yield at normal incidence is less than one, but 
for glancing incidence it is possible for values of N pos = 2 or 3 to occur (see bottom 

figure on page 14). Then from equation (2.35) we find 

E s max 2 = Ep ' E sl ' § (2.37) 

and so on, and proceed in the same way as for insulators. 


D. Reflection of Primary Electrons at Low Energies 

When a primary electron enters an insulator, the most probable type of collision 
is elastic scattering. When E p < E 0 , the threshold energy, the only kind of scattering 
process that can occur is elastic scattering. Even when E p is somewhat greater than 
Eq, inelastic scattering events have smaller probabilities and elastic scattering still 
tends to dominate. The Monte Carlo simulations for low E p primaries often show a 
track in which the primary electron goes through anywhere from a few up to more 


than a hundred collisions before eventually escaping from the surface. In these 
cases the energy of the escaping electron must still equal E p , and we speak of elastic 

backscattering or reflected primary electrons. 


Fig. 1 1 on the following page shows the backscattering yield, or number of 
elastically reflected primary electrons per total number of primary electrons, for 

E ranging from 0 to 50 eV. It will be seen from the figure that backscattering is quite 

important for low energy primaries from 0 to 15 eV. On the figure, the X's are the back- 
scattering yield obtained from the Monte Carlo simulations for argon. To this data we 
have fitted the following empirical formula: 


5,(E P < 16.8) = A E p b e- (E V c)d 
5 r (E p > 16.8) = 0.06 


The values of the coefficients as determined by least squares fitting are: 

A = 0.40 c = 8.1 

b = 0.70 d = 2.0 

This form gives a rapid rise from E p = 0, a maximum less than one at E p * 5.0 eV, 
a rapid decline toward E p = 16 eV, and then a constant value for Ep > 17 eV. The 
fiqure shows the empirical curve of equation (2.38), and the O’s show the result of a 
trial run with the average yield at each energy determined by using random numbers. 
The results are in good agreement with the general behavior of Rezvani s simulations. 

Since a reflected primary has undergone a number of collisions inside the solid 
before escaping from the surface, it has lost track of its original angle of incidence. 

(See Fig. 4.) Thus this phenomenon differs markedly from specular reflection, in which 
the angle of reflection is equal to the angle of incidence. We have looked at the angu- 
lar distribution of reflected primaries at several energies, and find that the results for re- 
flected primaries are similar to those for “true” secondaries. Therefore we can use the 
same angular distribution described in the next section for reflected primaries as well 
as for true secondaries. 

When the plasma-insulator simulation tells us that an electron strikes the sur- 
face of the insulator, we proceed as follows. First we calculate the reflection yield 5 re) 
from the empirical formula above, and then determine an integer yield in the usual 
manner. 


N« n = iNT(8 ren + r) 


(2.39) 
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Fig. 11 
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E. Angular Distribution of Secondary Electrons 


The task of assigning directions in spherical-polar coordinates is somewhat 
simpler than assigning energies to the secondaries. By examining Rezvani’s simula- 
tion results, we find that there exists azimuthal symmetry about the point of emission, 
as expected physically. Thus the azimuthal angled can be chosen randomly between 
0 and 2k. 

$ = 271 X Ot (2.40) 

where a is a random number such that 0 < a < 1 . 

Next we need to select the polar angle 9, measured from the outward normal to 
the surface (see Fig. 13 on page 31). For this purpose we plot results from Rezvani’s 
simulations on solid argon, as shown in Fig. 12 below. For these plots it is convenient 
to introduce a new variable x = cos 0. Each distribution is normalized properly by 
using the basic definition of a distribution function, namely, 

AN = N f(x) Ax (2 . 41) 



Fig. 12 Angular distribution of secondary electrons. The data points are 
simulation results for SEE from solid argon. Dot-dash line and O’s 
are for E p = 400 eV. Dash line and X’s are for E p = 50eV. The 
solid line is the best fit to the combined sets of data. 
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where N is the trial number ot events, and AN is the number of events in the range of 
width Ax. Then the normalized probability is given by 



1 AN 
N AX 


(2.42) 


There is a moderate amount of scatter in the simulation results because Simula- 
tions were mrf for only 100 primary electrons at each value of E p . But within these l,m- 
its there is no significant difference between the angular distributions for the two very 
Smeren. values of E p a. 50 eV and at 400 eV. Therefore we combine the resu ts o 

both distributions for the purpose of fitting an angular distribution function >" order 
a^ldate some of the culture, we perform a least squares fit with a function of the 

form: 


f(0) = A cos 0 + B cos 2 0 + C cos 3 0 


(2.43) 


Since only those electrons that make an acute angle with the outward normal 
will be emitted (Fig. 13), we normalize the distribution for 0 ranging from 0 to rc/2: 

/• tc /2 


rn/ 2 


1 = P(|) = 


dp(0) = N 


o 


f(0) sin(0) d0 


(2.44) 


= N 


rn/2 


0 


(A cos(6) + B cos 2 (0) + C cos 3 (0)} sin(0) d0 


(2.45) 


With the substitution x = cos 0 we obtain the result 


1 = N 


r i 


Jo 


(Ax + Bx 2 + Cx 3 } = |A + j + j 


N = 


I A + B ,CH 
12 3 }i 

H 


(2.46) 


(2.47) 
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We determined the coefficients A, B, and C by fitting this distribution function 
to the combined sets of simulation data for 50 eV and for 400 eV primary electrons. 
The result is that A = 0.9900, B = 1 .2375, and C = 0.0263. A plot of the best fit line is 
shown along with the data in Fig. 12. 

In most of the literature on SEE an angular distribution proportional to cos 0 is 
assumed. Then the normalized distribution function would be : 


f(0) = 2 cos 0 


(2.48) 


This is a simple form and is qualitatively justifiable on the basis of physical arguments. 
If we consider a cloud of internal secondaries, for instance, just beneath the surface, 
the situation is similar to a cloud of molecules in a gas near some surface. The 
velocity distribution of molecules actually striking the surface is weighted with a factor 


V n = V COS 0 


(2.49) 


where v n is the normal component of velocity. A similar distribution should be reason- 
ably appropriate for secondary electrons that manage to escape from the surface. 

But the detailed processes of secondary electron emission are more complicat- 
ed than this simple picture. The figure on page 6 of Part I shows the jagged and ir- 
regular paths followed by each secondary electron in the solid. Furthermore, the in- 
ternal secondary electrons are not in thermal equilibrium with the solid or with each 
other, as mentioned previously. There is even a potential barrier to overcome at the 
surface, so it is even less likely for a secondary electron with a small v n to escape. We 
believe that this is the reason that our angular distribution, as fitted to the Monte Carlo 
SEE results, is slightly more peaked in the forward direction than a simple cos 6 distri- 
bution. 


We have examined angular distributions for additional Monte Carlo simulations 
of SEE with Rezvani’s programs, especially for low primary energies of 10 to 50 eV. In 
all cases we find a good fit with the angular distribution function described here. Even 
for reflected primary electrons, which predominate for E p below 10 eV, we find the re- 
flected electrons follow a similar angular distribution. This result occurs because most 
of the reflected primaries make a number of elastic collisions (5 to 200) inside the 
solid before escaping from the surface. So we can use the same angular distribution 
function for all cases of reflected primaries and true secondaries. 
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OnrP the anaular distribution function is specified, we can develop the proce- 

Then, as for assigning energies to “| U9 of the integrated distri- 

,h3t “ITSS v integrated distribution, at equal intervals of integrated 
inverse look-up table nsiing u ■ « onnrnnriatp anale 0 bv table look-up and in- 

results shown at the end of this report. 



Fig. 13 A secondary electron being emitted from solid bulk material. 
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F. Determining Velocity Components of Each 
Secondary Electron 


When all the secondaries have been assigned an energy and direction, the 
only remaining task is to determine the components of velocity for each secondary 
electron. Recalling equation (2.7), we find the magnitude of the secondary electron’s 
velocity as: 


2 xE s x 1.6021 89xl0 19 \ I/2 , 

Vs = l BE 1 m /s 


(2.50) 


Then the velocity components are given by (Fig. 13): 


V x = V x sin(9) x cos(<(>) 
V y = V x sin(0) x sin(4>) 
V z = V x cos(0) 


(2.51) 
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III. Concluding Remarks 


The procedures outlined in Part II and the FORTRAN code m Appendix C were 
developed to realistically simulate secondary electron emission in bulk solids caused 
bv'incident electrons in a space environment. We have made no attempt to include 
secondary electron emission for materials in any other physical form, or emission in- 
duced by^other particles or photons, or high temperatures. The interested researcher 
is referred to the literature for such cases. Our model will produce realistic results 
primary electron energies between 0.0 eV and 1 .0 keV with an angle of incidence be- 
tween 0.0 and n/2 radians. However, some words of caution are in order. 

First the experimental and simulation data we used to develop our distribution 
functions and yield curves are generally from coarse data with the first data point at 50 
oMOO eV^ and the highest at 1 0 or 1 .2 keV. Experimental data for low energy primary 
electrons collected by Bruining show fine structure in the secondary electron yield 5, 
whereas the smooth curve obtained from equation (2.1) is actually an extrapolation 
back to 0.0 from 50 eV. But since the shape of the fine structure is dependent on the 
type of material and the fine structure in the actual yield curves produces little variation 
from equation (2.1), we decided to ignore such features in order to gain breadth of ap- 
plicability. By running some additional simulations with Rezvam’s programs at low E p 
(0 to 40 eV) we have been able to describe the behavior of SEE in this region more 
realistically. ’ The description of section IIC now gives a good picture of the phenomena 
involved with low energy primaries. 


Second, Bruning suggests that as the primary electron energy increases from 
0.0 eV to some small value the probability that the primary electron suffer sspecular re- 
flection increases. But the discussion of this in the literature is sparse, and little quanti- 
tative information is available. Fortunately the new runs with Rezvam’s programs 
demonstrate the behavior of reflected primaries (elastically backscattered primary 
electrons) very nicely. The large yield of reflected primaries for Ep < 10 eV, as de- 
scribed in section IID, will have a noticeable effect in the plasma simulations for solar 
panels. 


Third, the secondary electron energy distribution function we developed effec- 
tively goes to 0.0 at 30.0 eV with 8-bit floating point numbers. Our distribution is cer- 
tainly good for simulations where the primary electron energy is less than 1.0 kev. 
However at higher energies the number of inelastically backscattered primary elec- 
trons increases, but our distribution (and Rezvani’s model) does not account for the 
characteristic elastically backscattered primaries at larger primary energies 
(E p >100 eV). If these cases are to be accounted for we suggest development of 
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another distribution that will account for the required higher secondary energies. Then 
the programmer can choose between our energy distribution or the new one in an IF 
block as needed. 

Fourth, the literature generally subscribes to a cosine distribution for the angular 
distribution of the secondary electrons. This form was used by early researchers in 
this field and has remained the traditional distribution. We, on the other hand, have fit 
the distribution function of eq. (2.22) to the simulation results. There is actually not a 
great deal of difference between our angular distribution and a cos 0 distribution. We 
believe that our distribution is somewhat more realistic in representing enhanced 
emission at forward angles, near the normal direction. 

The plots on the following pages are from several runs of our simulations at 
designated energies of the primary electron; E p is set in the routine SECONDARIES 
and 0 p is maintained at zero. All the simulations were run with a sum total of less 
than 35 seconds of CPU time. 

All the plots show the behavior for model quartz, with band gap E g = 8.0 eV. 
The first three figures give the yield 8 versus primary energy E p from 0to50eV. 

The influence of the finite band gap in reducing the yield of true secondaries for 
Ep < 15 eV is evident in the first plot. The second plot shows the importance of reflect- 
ed primaries for E p < 15 eV, and the third plot gives the combined results for all eleo 

trons coming back from the surface. 

The next four figures show the energy distribution of the secondary electrons. 

At low primary energies (E p = 9 and 14eV) the influence of the band gap on the 
energy distribution is readily apparent. At larger primary energies (30 and 50 eV), 
the only noticeable effect is the backscattered peak at E s = E p . 

The angular distribution of the secondary electrons is shown in the last two 
figures. The distribution for the polar angle is peaked somewhat toward the normal 
(9 _ o°, cos 0 = 1), as discussed previously.. The distribution for the azimuthal 
angle is flat, as expected, with p(<J>) = p 0 = 1/(2 tc) = 0.159. 

As can be seen on the plots, the simulation results follow well the analytic distri- 
butions and arguments of Part II. Providing correct results with a minimum of CPU time 
is the best that can be expected of a numerical simulation of this nature. For this we 
have strived and, within the limits of our model, this we have obtained. 
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Fig. 14 Yield of secondary electrons vs. primary electron energy. The X’s are from Rezvani's simulation 
of the SEE process for quartz. The solid line is from the Haffner formula fitted to these results. 

The O’s show results from our procedures as described in the text, where the yield of secondary 
electrons at low primary energies is limited by the finite band gap of the insulator. 
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Fia 1 5 Yield of reflected primary electrons vs. primary electron energy for quartz. The symb 
9 ‘ and procedures are described in Fig. 14. Here the solid line is from the empirical formula 

for reflected electrons as given in the text. 
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Fig. 16 Total yield of secondary electrons for quartz. These are the combined results for secondary electron 
emission from Fig. 14 and for reflected primary electrons from Fig. 15. The O’s show the results 
from our procedures, and the solid line is from the Haffner formula. 
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RESULTS FROM 1000 PRIMARIES AT 14.0 eV 



Fig. 17 Energy distribution of secondary electrons 

for primary energies of 9.0 eV and 14.0 eV. 
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RESULTS FROM 1000 PRIMARIES AT 50.0 eV 



Fig. 18 Energy distribution of secondary electrons 

for primary energies of 30.0 eV and 50.0 eV. 
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Fig. 19 Angular distribution of secondary electrons. The upper distnbution fo 

the polar angle 0 is slightly more peaked in the forward direction than 
a simple cos 6 distribution. The lower distribution for the azimuthal 
angle <}> is uniform as expected. 
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C. FORTRAN Programs 


c subroutine SECDATA 

c 

c This subroutine is used to input the data files for secondary electron 
c emission. It is used only once at the beginning of a simulation. 

c sec_energy_file2 the data of secondary energy vs. probability 

c secjheta Jile the data of secondary angle vs. probability 

c sec_energy_file the data of probability vs. secondary energy 

c 

subroutine secdata 

common /energy/ AREAeng(3001),VALUEeng(3001) 
common /angle/ AREAthe(501),VALUEthe(501) 
common /energy2/ AREAeng2(1001),VALUEeng2(l001) 
c 

open(unit=92,file='sec_energyJile2.dat',status='old') 

open(unit=9l ,file='sec_theta Jile.daf, status=’old') . 

open(unit=90,file='sec_energy_file.dat',status='old‘) 

do i=1 ,3001 

read(90,*)VALUEeng(i),AREAeng(i) 
if(i.lt.502)then read(91 ,*)AREAthe(i),VALUEthe(i) 
end if 

if(i.lt.1002)then read(92,‘)AREAeng2(i),VALUEeng2(i) 

end if 

end do 
close(unit=90) 
close(unit=91) 
close(unit=92) 
c 

return 

end 

c subroutine SECONDARIES 

c 

c Secondary Electron Emission 

Q 

c Given on input the primary electron's velocity components [Vx], [Vy], and [Vz], and 

c an integer value [MAT] to indicate the type of material involved, and an integer 
c value [BOUND] to indicate the boundary surface, this routine will calculate the 

c primary electron's kinetic energy [Ep] and its angle of incidence [THETAP]. This 
c routine will also determine the probability of the primary electron elastically 
c reflecting with a call to "YIELDMAX". 
c 
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c Afterwards, the secondary electron production calculations will begin with a call to 
c "YIELDMAX," which will calculate a secondary electron [YIELD] based on the 
c Haffner formula. This real number is added to a random number less than one and 
c the sum is truncated to an integer. This integer is considered to be the total 
c possible secondary yield [Npos]. 

c This routine will then call a second routine called "SECENERGY" to assign an 

c energy [Es] to each secondary electron. To conserve the energy delivered to the 
c material by the primary electron, [Ep] + Ki will be considered the maximum amount 
c of energy available to all the possible secondaries. Each ( Jth) secondary will 
c receive an energy from a Monte Carlo method until there is no longer enough 
c energy to allow the (Jth+1 ) secondary to overcome the material's surface barrier, 
c or until J = [Npos]. 

c Then a third routine called "SECANGLES" will be called to assign the angles 
c [THETAS] and [PHES] to each secondary that had enough energy to overcome 
c the surface barrier. [THETAS] will be choosen with a Monte Carlo method first, 
c then [PHES] will be choosen at random, 
c 

c With the angles and energies of each secondary now assigned, this routine will call 
c its final subroutine "SECVELOCITIES" to assign velocity components [Ux], [Uy], 
c and [Uz] that are consistent with each secondary's assigned angles and energy, 
c 

c With all calculations now complete this routine will return the values [Nact], which 
c is the adual secondary electron yield, and the velocity components [Ux], [Uy], and 
c [Uz], which are arrays of dimension [Nad], 
c 

c If [Nad] = 0 then no secondaries were produced, but the primary stuck to the 
c material. If [Nad] = -1 then no secondaries were produced, and the primary was 
c refleded from the material. A postive value for [Nad] indicates the adual number 
c of secondaries produced, and that the primary stuck to the material. 

c 

subroutine secondaries(l,Vx,Vy,Vz, mat, bound, Uxx.Uyy.Uzz, Nad) 


c 


c 


c 


real Vx,Vy,Vz,thetap,Ep,Ki 
integer bound, mat 
integer I 


! Values of the primary eledron 
! Indicates the type of material 
! The passed random number argument 


real Uxx(30),Uyy(30),Uzz(30) 
real thetas(30),phes(30) 
real Es(30) 
real yield 

integer Nad, Npos 


! Secondary velocity 
! Secondary angular diredton 
! Energy of each secondary 
! Value from the Haffner formula 
! Adual and possible yields 


c 


logical reflect 
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c assign values to a, b, c, and Eg for Haffner formula 
c 


if(mat.eq.1)then 

! The material is QUARTZ 

a =1.02 


b =5.47 

! Eg is the band gap energy 

c = 14.02 


Eg = 8.0 

! Wf is the work function 

Ki = 1.0 


Wf = 0.0 


else if(mat.eq.2)then 

! The material is 

a =2.0 


b = 15.0 

! Holmes-insulator 

c = 4.72 


Eg = 4.0 

! Ki is the electron affinity 

Ki= 1.0 


Wf = 0.0 


else if(mat.eq.3)then 

! The material is Holmes-metal 

a = 0.1 


b = 10.0 

! a, b, and c are parameters 

c = 2.12 


Eg = 0.0 

! for the Haffner formula 

Ki = 0.0 


Wf = 5.0 


else if(mat.eq.4)then 

! The material is ALUMINIUM 

a =0.5 


b =8.3 


* c = 1.23 


Eg = 0.0 


Ki = 0.0 

! This value of Wf is from 

Wf= 4.25 

! Rezvani thesis, Vol 1 , page i 


end if 
c 

c calculate the energy and angle of incidence of the primary 
c 

Vsqrd = «Vx ** 2) + (Vy ** 2) + (Vz ** 2)) 

energ = 0.5 * 9.1 09534E-31 * Vsqrd ! Ep in JOULES 

Ep = energ/(1 .6021892E-19) lEpineV 

c 

if ((bound .eq. 1) .or. (bound .eq. 3) .or. (bound .eq. 5)) then 
thetap = atan(sqrt((Vx ** 2) + (Vy ** 2))/-Vz) 
else if (bound .eq. 2) then 

thetap=atan(sqrt(Vx**2+Vz**2)/Vy) 
else if (bound eq. 4) then 

thetap=atan(sqrt(Vx*‘2+vz**2)/-Vy) 

end if 
c 

c write(16,7 ’ 

c write(l6,*)'ep=\ep,' thetap=',thetap 


c calculate the possibility of reflection 

C if (mat .EQ. 1 -OR. mat .EQ. 2)then 

call elasref(l,Ep, reflect) 
if(reflect)then 

Nad = 1 


! Reflection 

! The refledion is elastic. 
! but not specular 


Es(l) = Ep 

call secangles(l, Nad, thetas, phes) 

call secvelocities(bound, Nad, Es, thetas, phes, Uxx.Uyy.Uzz) 


Nad = -1 


return 


end if 


end if 


0 

c calculate the total possibile yield 
c 

y = ran(l) 

call yieldmax(a,b,c,Ep,thetap, yield) 

Npos = int(yield + y) 


! This is the integer value of 
! secondaries that may be produced. 


c assign an energy to each secondary and find [Nad] 

call secenergy(l,Ep,Eg,Ki,Wf,Npos,Es,Nact) 


c 

c 

c 

c 


The remaining energy me, be considered to have been absorbed by the cystal 
lattice. Thus Ep is now zero and the primary is STUCK. 


if(nACT.eq.O) then 
return 

end if 


j Even though there existed the possibility of 
! producing secondaries, because of energy 
! considerations none were produced. 


c assign the diredional angles to each secondary 

call secangles(l, Nad, thetas.phes) 

c calculate the components of the velocities for the secondaries 

call secvelocities(bound,Nact,Es, thetas, phes, Uxx.Uyy.Uzz) 

c all calculations are complete for that primary electron 

c t 

c write(16,*)'npos=',npos/ nact= ,nact 

c do i=1,nact 

c write(l6,*) , i=',i.' es=',es(i),’ thetas' ,thetas(i) 

c end do 

return 

end , . 

c=========================e nd °* secondaries========== 


c subroutine SECANGLES 

c This routine will assign values to the spherical-polar angular 
c coordinates theta and phe for each secondary eledron. 


c The azimuthal angle, theta, will be generated with a Monte 
c Carlo method using the distribution function: 
c 

c Y = A COS(theta) + B COS A 2 (theta) + C COS A 3(theta) 
c 

c where the parameters A,B, and C are determined by a best fit to G. A. Rezvani's 
c secondary electron angular distribution data for runs at 50 eV and 400 eV. 
c The polar angle, phe, will be chosen at random since the determination of 
c phe = 0.0 is arbitrary to the plasma simulation. 

c=== ================================================== ============= 

c 

subroutine secangles(l,Nact, thetas, phes) 
c 

integer Nact 

real thetas(30),phes(30),rtheta 
c 

c find the angles 
c 

do i=1,Nact 

rtheta = ran(l) 

call find_theta(rtheta,thetas(i)) 
phes(i) = ran(l) * (2 * 3.141593) 

end do 
c 

return 

end 

c subroutine SECVELOCITIES 

c 

c This routine will assign values to the velocity components 
c [Ux], [Uy], and [Uz] for each secondary electron, [Vs], 
c 

c For each secondary electron the magnitude of the total velocity [V] is calculated 
c from a knowledge of its total energy [Es]. Then each component is calculated 
c according to the following equations: 
c Uz = V * cos(theta) 

c Uy = V * sin(theta) * sin(phe) 

c Ux = V * sin(theta) * cos(phe) 

c 

subroutine secvelocities(bound, Nact, Es, thetas, phes, Uxx.Uyy.Uzz) 
c 

integer Nact.bound 

real Es(30),Uxx(30),Uyy(30),Uzz(30),thetas(30),phes(30) 
real m,V 
c 

m = 9.109534E-31 
c = 1.602189E-19 


do i=l ,Nact 


V = sqrt((2 * Es(i) * c)/m) 

Uxo = V * sin(thetas(i)) * cos(phes(i)) 

Uyo = V * sin(thetas(ij) * sin(phes(i)) 

Uzo = V * cos(thetas(i)) 
c 

if ((bound.eq.1).or.(bound.eq.3).or.(bound.eq.5)) then 
Uxx(i)=Uxo 
Uyy(i)=Uyo 
Uzz(i)=Uzo 

else if (bound. eq. 2) then 
Uxx(i)=Uxo 
Uyy(i)=-Uzo 
Uzz(i)=Uyo 

else if (bound. eq.4) then 
Uxx(i)=Uyo 
Uyy(i)=Uzo 
Uzz(i)=Uxo 

end if 

end do 

return 

end 


c========================= end of secvelocities 

c================= ======================== 

c subroutine SECENERGY 

c 


c This routine assigns an energy [Es] to the secondary electrons. The information 
c that comes from the calling routine is the total available energy [Eavail], the 
c threshold energy from the material's band gap [Eg] and electron affinity [Ki], and 
c the total possible yield from the (modified) Haffner formula. This routine will loop 
c until all the electrons given by the Haffner formula have been assigned an energy, 
c or until the energy equation cannot be satisfied: 
c 

c Eavail = > Ki + Eg + Es(l) + sum[Es(l-1 ) + Eg + Ki] (1 ) 

c 

c Since each secondary must overcome the surface barrier, it will lose an amount of 
c energy equal to Eg + Ki before it is outside the material. The sum of the energies 
c assigned to each secondary plus the sum of the energy lost as each secondary 
c crosses the surface barrier must not be greater than the energy delivered to the 
c material by the primary electron. That way energy is conserved. After going 
c through this loop as many times as possible any "left over" energy will be 
c considered to be absorbed by the crystal lattice, again conserving energy, 
c 

c The distribution function that will be used to assign secondary energies is from 
c the best fit to G. A. Rezvani's data on secondary energy: 
c 
c 
c 


f(Es) = A(Es A 0.5)exp[-(Es/Eo) A 2] (2) 


c A random number will be chosen and will be associated with the integral of the 
c distribution function from zero to some Es, which is normalized for energies 
c between 0 and some [Emax]. Once this area is known, a call to "FIND_ENERGY H 
c will give the value of Es(l). If Es(l) doesn't satisfy equation (1), Ep is considered 
c to be exhausted. 

c 

subroutine secenergy(l,Eavail,Eg,Ki,Wf,Npos,Es,Nact) 
c 

real Eavail, Eg, Ki, Wf, Eo, Es(30) 
realdlEmax, Emax, rarea 
integer Npos, Nact 
c 

c calculate Es and the actual_yield 
c 

Nact = 0 

Eo = Eg + Ki + Wf 
c 

dIEmax = Eavail + Ki + Wf 

if (npos .gt. 30) npos=30 
do i=1 ,npos 

Emax = dIEmax - Eo 

if(Emax.lt.O)return 

rarea = ran(l) 
call find_energy( rarea, 
c 

dIEmax = Emax - Es(l) 

Nact = i 

end do 
return 
end 

c=========================end of secenergy=========================== 

c subroutine FIND_ENERGY 

c 

c This routine selects an energy for the secondary in accordance with 
c the distribution function described in the text. 

c 

subroutine find_energy(proba, Emax, eng) . 
c 

real proba,Emax,eng 

common /energy/AREAeng(3001),VALUEeng(3001) 
common /energy2/ AREAeng2(1001),VALUEeng2(1001) 


! The magnitude of the surface barrier. For 
! insulators, Eo is the sum Eg + Ki. For metals, 
! Eo is just the work function. 

! The maximum energy deliverable to 
! interactions with the secondaries. 


I The total possible energy the secondary 
! under consideration can have. 

! There is not enough energy to get 
! another secondary past the surface barrier. 

Emax, Es(i)) 

! The total deliverable energy less the amount 
! delivered to the secondary being considered. 
! An actual secondary has been observed. 


probal =proba 
if (Emax .It. 30.0) then 

j=int(Emax*100)+1 

Pmax=AREAeng(j)+(Emax-VALUEeng(j))*(AREAengG+1) 

& -AREAeng(j))/(VALUEeng(j+1)-VALUEeng(j)) 

probal =proba1 * Pmax 

end if 
c 

index=int(proba1*1000)+1 

eng=VALUEeng2(index)+(proba1-AREAeng2(index))‘(VALUEeng2(index+1) 

& -VALUEeng2(index))/(AREAeng2(index+1 )-AREAeng2(index)) 

c 

return 

end 

-____ ====== ===============end of find_energy=============== = - 


c subroutine FIND_THETA 

c 

c This routine selects a value of theta for the secondary electron, 
c 

subroutine find_theta(proba, angler) 
common /angle/ AREAthe(501),VALUEthe(501) 
c 

index=int(500.0‘proba)+1 

angler=VALUEthe(index)+(proba-AREAthe(index))*(VALUEthe(mdex+1) 

& -VALUEthe(index))/(AREAthe(index+1)-AREAthe(index)) 

return 

end f d h t 

c subroutine YIELDMAX 

c 

c This program will make use of the angular modified Haff ner formula: 

c Yield(E) = C[EXP(-aE) - EXP(-bE)]*EXP(1 - COS(THETA)) 

c to calculate the secondary electron [YIELD] as a function of the primary electron 
c [ENERGY] and angle of incidence [THETA]. 

c 

subroutine yieldmax(a,b,c, energy, theta, yield) 
real yield, energy, c, a, b, theta 
c 

energy = energy/1 000 ! convert Ep from eV’s to keV’s 

aterm = exp(i-cos(theta)) 

yield = c * ((exp(-a * energy))-(exp(-b ‘energy))) * aterm 

energy = energy * 1000 ! convert Ep back to eVs 


s 


return 

end 


===end of yieldmax 


c subroutine ELASREF 

c This routine will take an empirical formula fit to G. A. Rezvani's ; simulation idata i for 
c elastic reflection of the primary electron at a broad number of energies tetwee 

c 1 0 and 50 0 eV. Since the result of the equation is always a fraction and since the 

c primary must reflect or not reflect on an all or nothing basis, we will add a random 
c fraction to the result of the lormula, take the integer part and truncate it to an 
c integer If the result of this proceedure is one, we consider the primary electron to 
c have suffered elastic reflection. If the result is zero, there is no reflection. 

c 


subroutine elasref(l,Ep, reflect) 
real Ep 

integer l.i.ichance 


logical reflect 
reflect = .false. 

a = 0.40 
B = 0.70 
C = 8.1 
d = 2.0 


! These are from a fit to Rezvani's data 
! for Ep equal to or less than 1 7 eV. 


if(Ep.lt.16.88)then 

y = a* (eP*‘b) * exp(-(eP/c)“d) 

else 

y _ o 06 I At higher energies, Rezvani's data flatten out 

end if 

ichance = int(y + ran(l)) I The chance to reflect 

if(ichance.gt.0)reflect = .true. 


return 

end 


======end of elasref======== 
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